Dynamic analysis of lathe bed of woodworking CNC machining center based on the modeling of joint surface

Based on the finite element theory, a joint-plane modeling method is employed to connect the corresponding nodes at the joint surface of the woodworking computer numerical control (CNC) machining center bed with a 2-node 12-degree-of-freedom unit. A spatial element model is established, which can show the state of the nodes between joint surfaces when they are stretched, compressed, or twisted; and it can help build a woodworking CNC machining center on a finite element model of bed with the characteristics of the joint surface. The simulated analysis is performed on the model and is compared with the result of simulated analysis on the bed model that ignores the characteristics of the joint surface and modal experiment. The comparison verifies the effectiveness of the modeling method based on the characteristics of the joint surface. The weak link of the machine bed structure is analyzed and optimized. The natural frequency of the bed is improved by2.55% ~ 11.3%. The displacement is reduced by a maximum of 19.4%, and dynamic performance of the bed is improved.


Introduction
Modern industry can provide a highly flexible, individualized, and digital production and service. With the standardized protocol, interoperability between networked machines is realized and dynamic reorganization of machinery and modules in the factory is achieved [1]. Woodworking, the same as other industries, is entering the era of Industry 4.0, conducting digital production and processing. To meet the personalized needs of the product, the woodworking deeply integrates information technology with manufacturing technology, forming a woodprocessing intelligent system based on the industrial Internet of Things (IoT). Therefore, intelligent manufacturing is an inevitable trend in development of woodworking machinery. The processing technology and production function are continuedly improved alongside with profound integration of information and manufacturing technologies, giving higher demand for dynamic performance of woodworking machinery to ensure maximum productivity, output, and reliability of the overall manufacturing process [2]. Woodworking CNC machining center, as the core equipment in IoT of the digital wood production and processing, can connect other mechanical equipment to achieve an equilibrium quantity and optimize cost in the possibly shortest time to efficiently produce highly personalized products. Its dynamic performance influences the efficiency and reliability of the entire intelligent production system.
In the mechanical structure, the surface in which the parts and components contact each other is a mechanical joint surface. Performance of the joint surface can affect the dynamic characteristic of mechanical structure. 30%-50% of the machine tool stiffness depends on the stiffness of the joint surface; more than 90% of the machine tool damping value is derived from the damping of the joint surface; and 60% of the machine vibration originates from the characteristics of the joint surface [3]. Characteristics of the joint surface can be affected by many factors, most of which are nonlinear factors. Various processing methods and application conditions are intertwined, resulting in very complex characteristics of the interface. For the complexity, establishing a theoretical model that accurately characterizes the stiffness and damping of the joint surface is the key to establishing a finite element model of the mechanical structure, and it is the premise for predicting and analyzing the dynamic characteristics of the mechanical structure.
Some specialists and academics have conducted extensive research on the subject of the dynamic characteristics of machine tools. Hao Y. P. et al. proposed a multi-source signal fusion flutter detection approach based on WPD and power entropy, which can more accurately detect the existence of early flutter and predict its severity [4]. Zhu L. D. et al. analyzed the research progress of milling flutter prediction, detection, and suppression, and the research proposed that integrating flutter prediction, detection, and suppression units into intelligent machine tools or intelligent spindles is the future path of research [5]. In addition to studying the main vibration characteristics of machine tool structure under cutting conditions, Huang Q. et al. proposed a method to predict the modal response signal and identify the dominant vibration mode of machine tools [6]. The aforementioned research results indicate that the dynamic characteristics of machine tools have a significant effect on the precision and chatter of the machining system, which must be studied further to improve the quality and efficiency of machining.
Some experts and scholars mainly focus on the experimental research and analysis of the dynamic characteristics of the whole machine and various mechanical system structures of woodworking machinery. Pimenov D. Y. et al. [7] studied the dynamic characteristics of woodworking CNC machining center and analyzed the amplitude-frequency characteristic of the CNC machining center and the dynamic performance of the woodworking CNC machining center during external excitation and milling. Noda [8] used the bending theory to derive the equation of motion under the action of thermal load and plastic strain of the tensioned circular saw. To study the influence of tension conditions on the dynamic characteristic of tensioned circular saws and find out the appropriate tension condition; Marinenov [9,10] analyzed the deformation of the large band saw machine spindle and transmission system caused by the static and dynamic load; and parametric equations reflecting the deformation of each system are obtained to ensure safe and reliable operation of the machine related to various optimizations. The current research on dynamic characteristic and vibration of woodworking machinery generally ignores the influence of the characteristics of the joint surface, and few studies can combine theoretical analysis with experimental methods.
In the process of optimizing the machine tool structure, Shen L. et al. performed sensitivity analyses on each component based on the results of modal tests, chose the optimization target, and finished optimizing the spindle box, pillar, and lathe bed [11]. Based on the natural development law of leaf veins, Hao Y. P. et al. proposed a new method to design the architecture of reinforcement inside the supporting sections of machine tools [12]. The optimization objectives and design variables of the biomimetic CNC turntable were chosen by Lui S. H. et al. based on the results of thermodynamic coupling and modal analysis. They then used the sensitivity analysis method to investigate the influence law of the design variables on the optimization objectives, found the biomimetic CNC turntable's optimal solution through response surface optimization, and completed the optimization design [13]. The aforementioned research results demonstrate that the current mechanical structure optimization design primarily uses mathematical programming as its core, is based on a number of research methodologies, and increases optimization accuracy while maximizing a greater number of objectives.
This paper offers a novel mechanical structure modeling method based on finite element theory that takes into account the dynamic characteristics of the joint surface. In this method, the mechanical structures of each part are combined and connected one-to-one corresponding nodes on the joint surface using two-node space elements. This simulates the tensile, compressive, and torsional characteristics between the mechanical joint surfaces. This modeling technique is used to investigate the dynamic properties of the bed of the PTP160PLUS Woodworking CNC Machining Center (hence referred to as the bed), construct a finite element model of the bed, and perform a modal analysis on the model. With the use of the hammer method, the bed is carried through its modal test, and the results of the test are used to validate the modeling method. On the basis of the results of the modal analysis, the weak links of the bed structure are identified, and the bed structure is adjusted using the sensitivity analysis method in order to enhance the dynamic properties of the bed.
The rest of the paper is organized as follows. In Section 2, the bed finite element modeling is discussed. In Section 3, the dynamic analysis and modal test of bed are discussed. The bed structure is optimized in Section 4. Finally, the conclusion is given in Section 5.

Modeling method of joint surface
The modeling at the joint surface is based on its dynamic characteristics, which is discretized and applied to the finite element model of the integrated composite structure. Literature [14] pointed out that stiffness, damping, and surface pressure of the joint surface show a continuous power function relationship, but the difference is not large. In the elastic range, the stiffness, damping, and surface pressure of the joint surface can be simplified to linear calculations. The finite element method is adopted to divide the joint mesh. If the mesh of the joint surface is sufficiently thin, characteristics of the joint surface will be equivalent to the joint between the joint nodes [15,16]. Therefore, this paper uses the joint surface nodes to establish the joint surface equivalent model the joint is evenly divided into enough regular shape elements and nodes while ensuring one-to-one correspondence of nodes on the bonding surface and between each pair of bonding nodes [17]. A set of 12-degree-of-freedom space spring damping unit is adopted to simulate the property between the bonding faces. The reasons for using this unit to express the properties between binding surfaces are as follows: 1. The nodes at both ends of the cell contain 6 degrees of freedom, which is consistent with the interaction between the binding surface, so the unit can accurately express the binding surface properties.
2. The interactions between the binding surfaces are dispersed to the units between the nodes to transform the nonlinear problem into a linear problem to facilitate the numerical calculation.
3. Each unit has no spatial position relationship and quality characteristics, and will not add interference factors to the whole system.

The interaction relationship between the units can be realized by the boundary coordination equation and the constraint equation between the nodes.
In the joint plane in any actual woodworking machine structure, there are generally six forms of forces at the joint surface [9] (Fig 1a). These forces include the positive force in the μ direction F μ , the shear forces F ν and F ω in ν and ω directions, the bending moments T ν and T ω around the ν and ω directions, and the torque T μ around the μ axis.
The finite element method is adopted to evenly divide the contact joint of an actual machine tool structure (Fig 1b), and obtain n sets of joint nodes. Any joint node i is selected on a certain structure of the joint surface, and there will be corresponding joint node j on the corresponding joint surface. A joint space stiffness damping unit is similar to the beam unit established between the joint nodes, and an equation at the joint surface is obtained according to the equilibrium equation of the joint surface [18]: where, F μ , F ν , and F ω are the loads in the μ, ν, and ω directions of the joint surface, respectively; T μ , T ν , and T ω are the torque and bending moment of the joint surface in the μ, ν, and ω directions, respectively; K μi and K θμi are the normal stiffness and torsional stiffness between the joint nodes; K νi and K ωi are the tangential stiffness between the joint nodes; K θνi and K θωi are the tangential torsional stiffness between the joint nodes; Δμ i is the normal deformation amount between the joint nodes; Δν i and Δω i are the amount of tangential deformation between the joint nodes; Δω μi , Δθ νi , and Δθ ωi are the torsional angular displacement between the joint nodes; and n is the number of joint nodes. For any group of spring damping elements in space that satisfy the tension, compression, and torsion characteristics between nodes in the joint region (Fig 2), the nodes that make up the element have 6 degrees of freedom to ensure the properties of the joint region, and the number of nodes that make up the space unit of this area is no more than 2. In the process of interaction, there will be elastic or plastic deformation between nodes [19], so the joint node satisfies the deformation relationship of two nodes, namely: where, [K] is the stiffness matrix between the joint nodes; [C] is the damping matrix between the joint nodes; F μi , F νi , and F ωi are the loads in the μ, ν, ω directions of the node i, respectively; F μj , F νj , and F ωj are the loads in the μ, ν, and ω directions of the node j, respectively; μ i , ν i , and ω i are the displacements of the joint node i in the μ, ν, and ω directions, respectively; μ j ,ν j , and ω j are the displacements of the joint node j at μ, ν, and ω directions, respectively; T μi , T νi , and T ωi are the bending moments and torques in the μ, ν, and ω directions of the joint node i, respectively; T μj , T νj , and T ωj are the bending moments and torques in the μ, ν, and ω directions of the joint node j, respectively; θ μi , θ νi , and θ ωi are torsional angular displacements in the μ, ν, and ω directions received by the joint node i, respectively; and θ μj , θ νj , and θ ωj are torsional angular displacements in the μ, ν, and ω directions received by the joint node j, respectively. When there is a joint structure modeling, the nodes of the joints have correspondence. Therefore, the stiffness matrix [K] and the damping matrix [C] of the joint unit are symmetric matrices. C ¼ where, There is no mass at the joint surface. According to the principle of continuous system discretization, the multi-degree-of-freedom system at the joint surface can be expressed as follows: where, {q} = [μ ν ω], K and C are the stiffness matrix and damping matrix of the joint, respectively, and F is the external load of the joint node. The positive direction of the vector is specified, and relationships between the stiffness matrix of each joint unit, the damping matrix, and the stiffness damping of the joint region are obtained from Eqs (2) and (3), respectively; and the model at the joint surface is established according to Eq (1). For any substructure: the finite element model can be established based on the dynamic equation (Eq (4)) of the continuous system of the structure itself. The model at the joint surface is coupled with the model of each substructure through common nodes and units to establish a machine tool finite element model with the characteristics of the joint surface.
where, {q r } = [μ r ν r ω r ]; M r is the mass matrix of the substructure; K r is the stiffness matrix of the substructure; C r is the damping matrix of the substructure; and F r is the external load of the node at the substructure. In addition to ensuring uniform joint surface of the bed and the correspondence between the joints, the accurate rigidity and damping of the joint surface of the machine tool can't be ignored in the overall modeling process of the bed. At present, the stiffness and damping of the joint are mainly obtained by identifying the test data and the model. The identified joint stiffness and damping are substituted into the stiffness matrix and damping matrix of Eq (2) to establish the joint unit model. The subunits are associated with each substructure model to establish a bed model containing the dynamic characteristics of the joint surface.

Establishment of finite element model of woodworking CNC machining center bed 2.2.1 Bed finite element analysis process.
Typically, the finite element analysis process consists of three phases. The first part is the model's pre-processing. The structure to be analyzed's geometric model, or 3D modeling, is established first. In this work, SolidWorks is used to create a 3D model of the bed, which is then meshes using a mesh division program. The second step is to specify the solution condition, which is the crux of finite element analysis and includes model selection, material creation, setting of constraint and boundary conditions, etc. In this study, HyperMesh is used to divide the mesh, Ansys is used to set the solution condition, and finite element solution is carried. The final part is post-processing. The post-processor can perform graphical display and statistics on the data obtained after solving in order to examine the calculation results, build and display the system displacement cloud map, and so forth. Fig 3 depicts the process of bed finite element analysis.

Bed 3D modeling.
When the finite element model is established, the 3D model of the bed is firstly established based on its actual structure. If the fine features such as rounded corners, chamfers, bolt holes, and small bosses in the bed are all reflected in the model, the number of cell grids is huge, the calculation accuracy is degraded, and the error of the simulation analysis is increased, greatly affecting the simulation calculation result [20]. Therefore, the bed structure can be simplified to improve the accuracy and accelerate the speed of the bed meshing and finite element analysis. Fig 4 is a simplified 3D solid model of the bed.

Bed model meshing.
In this paper, Hypermesh software is employed to mesh the bed model, and enough units and nodes with regular shapes are divided at each joint surface to ensure one-to-one correspondence of the nodes on the joint surface. In this way, it can establish the equivalent dynamic model of the joint surface through these nodes, and set the material properties of the various parts of the bed in the Hypermesh software, the material properties of each part of the bed are shown in Table 1, 458,231 units and 151,046 nodes are divided into the bed.

Joint surface modeling.
The main joint surfaces in the bed are between the bed and the support, between the bed and the slider, and between the support and the slider. There is a relationship between the stiffness at the joint surface, the damping, and the surface pressure of the joint surface. In the continuous power function relationship, according to the influence cone theory [21,22], the main concentration range of the surface pressure generated by the bolt pre-tightening force is the cone area near the circular gasket, and the angle of action is generally conservatively selected a = 30˚. Since the bolts on the bed are evenly distributed and dense enough, the stiffness and damping at the joint surface are more evenly distributed. According to Yoshimura's unit area integral method [23], the equivalent stiffness of the fixed joint surface can be calculated according to Eqs (5) and (6).
where, K 1 and K 2 are the normal stiffness and tangential stiffness of the joint surface, respectively; and k 1 (p) and k 2 (p) are the normal stiffness and tangential stiffness per unit area, respectively. The base of the machine bed and the ground are connected to the ground by bolts, and the simulation of the boundary conditions is defined by the fixed joint surface. The data per unit area of the fixed joint surface used in this paper is the basic data of the test. According to the bonding conditions of the joint surface (material, surface pressure, joint surface medium, roughness, etc.), the characteristic parameters of the unit area joint surface are selected and substituted into the Eqs (5) and (6) to calculate the equivalent characteristic of the actual joint surface ( Table 2). The surface pressure refers to the positive pressure per unit area of the joint  surface, which is calculated by the combination of bolt preload, component gravity, and cutting force. Similarly, the equivalent of damping can be completed. In this paper, the matrix unit Matrix27 is undertaken as the connection unit to integrate the characteristics of the joint surface. The Matrix27 unit is created on the joint surface generated during meshing. According to the parallel principle of the spring and the number of joint units [24], the characteristics of the joint surface are discretized to each unit, and the normal and tangential stiffnesses of the single unit on different joints are obtained, the command stream that used to add Matrix27 unit is shown in S1 Text. The Matrix27 element rigid and damp real constants are set to give the joint surface stiffness and damping property. The bed finite element model and joint space unit diagrams are shown in Fig 5.

Finite element modal analysis
The finite element modal of the bed-limited element model is analyzed by ANSYS. According to the actual working conditions of the woodworking CNC machining center, the modal is analyzed for the bed model considering the dynamic characteristics of the joint surface and the machine bed model with rigid connection of the joint surface. The block Lanczos method (Block Lanczos) is employed to extract the natural frequencies and modes [25]. Generally, the actual working speed and the frequency of the woodworking CNC machining center are5, 000 15,000 r/min and 83 Hz~250 Hz, respectively. Therefore, the first 6 th order modal data are taken according to the modal calculation result, and its natural frequency and mode shape description are shown in Table 3. f is the frequency of each step of the bed obtained by the modal test; and f 1 and f 2 refer to the analysis frequency obtained by the finite element simulation of the bed with and without the characteristics of the joint surface, respectively.

Test equipment.
To compare with the analysis frequency and vibration mode obtained by the limit element simulation, the validity of the characteristics of the joint surface modeling method is verified, and the bed structure is tested by the hammer method [26]. The test equipment mainly includes the pulse hammer, the acceleration sensor, the BK multi-channel data acquisition front end, and the ME'scope data processing system. The test modal test site is shown in Fig 6(a).

Test method and test point.
The single-point excitation multi-point response is adopted in the test. According to the result of finite element analysis, the point in the upper right corner of the workbench is selected as the excitation point. Based on the test point, structure of the whole system is reflected as much as possible without missing the principle of the mode [27]. 63 measuring points are arranged on the body. Fig 6(b) is a block diagram of the measuring point arrangement on the bed in the test. The intersection of each line is the position where the measuring point is.

Comparison of test and simulation modes and data
Frequency in the modal test is shown by f in Table 3, the error between the simulation and test result is shown in Fig 7, and the comparison between the mode vibration obtained by the modal test and that from the simulation analysis is shown in Fig 8. In Fig 7, the comparison of the simulation analysis and frequency in the modal test shows that resulting error of the bed model without characteristics of the joint surface is 13.28~3.08% larger than that of the bed model with, and latter is closer to reality. The first-order analysis frequency of the bed model with the characteristics of the joint surface shows a larger error of 19.55% compared with the modal test result, which may be because that the first-order vibration deformation mainly occurs on the slider structure. As the slider structure is too small, the test point can't be directly arranged on the slider structure during the modal test, so a nearby point is adopted instead. The error between analysis frequency of other order and test frequency is 1.08%~9.44%. According to Fig 8, the simulation analysis and the test order mode are almost the same. The comparison of conditions between the simulation analysis and the modal test suggests that the data error between the simulation and the test result may be caused by following two aspects: 1. The basic data of the joint surface in the modeling process is measured according to the test, and there is an error.
2. In the modal test, the energy input to the system by the hammer stroke is limited, so the accuracy of the signal received by the sensor decreases.
Despite of the numerical error between the simulation result and the test result, the data comparison suggests that in the bed modeling process, the bed model with the characteristics of the joint surface is closer to the modal test result, with a reasonable error. The above result shows that the bed model with the characteristics of the joint surface can better reflect the dynamic performance of the bed. This paper using the equivalent modeling method of the joint surface is proved effective, accurately reflecting the characteristics of the bed joint surface. In addition, the bed model established based on characteristics of the joint surface in this paper is proved reasonable and can be popularized and applied to other woodworking machinery for analysis of dynamic performance.
According to the analysis and test vibration mode result, vibration deformation of the bed mainly occurs on the workbench. The workbench can greatly affect the machining precision of the woodworking CNC machining center. It should be considered in design optimization and increase the thickness or ribs inside the workbench. The vibration mode of the bed  suggests that the front wall shows insufficient bending and torsion resistance, and the torsion and flexural rigidity can be improved by appropriately increasing the thickness of the sheet metal or arranging the ribs inside the bed.

Optimization of bed structure for the woodworking CNC machining center
According to the simulation analysis and experimental modal result, vibrations of the bed and the front wall are more obvious. The modal analysis is performed to optimize the bed structure, and the result of the finite element analysis of the bed model before and after optimization are compared to verify the optimization accuracy.

Establishment of the bed optimization model
According to the simulation analysis and experimental modal result, the weak part of the bed is optimized in the SolidWorks software. Two ribs are added to the inner wall of each workbench, two rib plates are added to the left and right sides on the front metal plate inside the bed, respectively. The optimized bed 3D solid model model is shown in Fig 9. In this optimization, considering the computation amount, only seven parameters are chosen as the input quantity: the distances L 1 , L 2 , L 3 between the end points A, B, C of the left side rib plate and the O intersection point of the bed's sheet metal, the thickness D of the rib plate, and the width W and height H of the rib. On the body mass, natural frequency, and relative displacement, the effects of these seven parameters are studied. Table 4 displays the range of values for each parameter.
The mass of the bed, the maximum relative displacement, and the first six natural frequencies are set as the objective functions for the structure optimization of the bed, and an optimization model is established. Fig 10 depicts the bed structure optimization process.

Optimization of bed structure
The sensitivity analysis method, which analyzes the sensitivity of structural performance function changes to structural design parameter changes, is used to optimize the bed structure. Using the sensitivity analysis function of the finite element analysis program ANSYS Workbench, the optimal parameter size can be easily identified. Since the reb on the workbench and the reb plate on the bed body act on distinct components, they can be taken into account individually during the optimization process to establish the optimal size. Fig 11 displays the sensitivity of the four reb plate parameters L 1 , L 2 , L 3 and D to the first six natural frequencies, as well as the maximum relative displacement and bed mass.
The software ANSYS Workbench is utilized to optimize the bed's structure, with the maximum sixth-order natural frequency, the least maximum relative displacement, and the maximum bed mass serving as optimization objectives. The solution is calculated using the MOGA genetic algorithm, and after multiple iterations, the optimization feasibility solution set is achieved, the solution set is listed in S1 Data. After thoroughly taking into account the four optimization objectives of the bed, three groups of relative optimal solutions Candidate A, Candidate B, and Candidate C are generated, as shown in Table 5.

PLOS ONE
The optimal choice for the best bed structure optimization solution is Candidate A following comparison analysis. The optimal solutions for the width W and height H of inner wall ribs in the workbench, using the same optimization solution process, are 7.94 mm and 10.11 mm, respectively. The parameters of each dimension are correctly rounded for the convenience of the real product design and manufacture, and the results are displayed in Table 6.

PLOS ONE
The optimized structural parameters are taken as the structural dimensions of all the rib plates and ribs of the bed, in order to establish the optimized finite element model of the bed, and the model is entered into ANSYS to do the finite element analysis.

Comparison of finite element analysis result before and after bed optimization
According to the bed optimization scheme, the reinforcing ribs are added to the inner wall of the workbench and two ribs are added to the inner wall of the bed near the front wall. Therefore, the total mass of the woodworking CNC machining center bed is only increased by 0.42%, and the bed structure changes not greatly before and after optimization. Comparison between the result of the meta-analysis and the natural frequency before the optimization is shown in Table 7.
Comparison on the result of finite element analysis before and after optimization of the bed structure reveals that the natural frequency of the optimized bed structure is increased by 2.55%~11.3%, and the corresponding relative displacement of each step is reduced with the maximum reduction value of 19.4%. It indicates that the optimized bed model shows increased rigidity, reduced vibration deformation, and greatly improved dynamic performance, achieving the design optimization.

Conclusion
1. The dynamic characteristics analysis of woodworking machine tools was conducted using a method based on the finite element method that took into analysis the dynamic characteristics of the joint surface. The equivalent model of the joint surface was established by directly simulating the tensile, compressive, and torsional characteristics between mechanical joint surfaces utilizing the nodes of the joint surface. Based on this method, the parameters for the joint surface were established in the finite element model of the bed of the woodworking CNC machining center.
2. The modal analysis results of the bed model with joint surface parameters were more compatible with the test results when compared to the simulation results of the model without joint surface parameters, combined with the modal test data, and the overall error was within 19.55%. It demonstrates that examining the dynamic characteristics of the joint surface is crucial for analyzing the dynamic characteristics of the bed, and that the modeling method of the joint surface presented in this study is effective.
3. The weak points of the bed structure were improved and optimized in accordance with the analysis of each mode form. The bed optimization model was established using the sensitivity analysis method, and the optimal size of the optimization parameters was determined. The natural frequency of the bed after optimization was enhanced by 2.55%-11.3% in comparison to the results of the simulation analysis before and after optimization, and the maximum related relative displacement was decreased by 19.4%. The dynamic performance of the bed was enhanced, and the intended optimization was accomplished.
Supporting information S1 Text. Command stream that used to add Matrix27 unit.